Two Stream Plasma Instability

Defining the Equations of Motion

In the two stream instability, we have a stream of ions $(j)$ stationary with respect to a stream of moving electrons $e$ in direction $x$ with velocity $V_0$. Thus, the equations of continuity become

\frac{\partial n_j}{\partial t} + \nabla \cdot n_j \mathbf{v} = \frac{\partial n_j}{\partial t} + \frac{\partial(n_j v_j)}{\partial x} = 0

\frac{\partial n_e}{\partial t} + \nabla \cdot n_e \mathbf{v} = \frac{\partial n_e}{\partial t} + \frac{\partial (n_e v_e)}{\partial x} = 0

where $n$ is the number density of the ion or electron stream, and $v$ is the velocity of the stream.

We assume that the plasma is both **cold** and **unmagnetized**, leading to the following momentum equations

n_j m_j \left[\frac{\partial v_j}{\partial t} + v_j \frac{\partial v_j}{\partial x}\right] = q_j n_j E

n_e m_e \left[\frac{\partial v_e}{\partial t} + v_e \frac{\partial v_e}{\partial x}\right] = q_e n_e E,

where $m$ is the mass of the individual particles and $q$ is the charge of the particle.

At this point we have four equations, but 5 unknown variables $(n_e, v_e, n_j, v_j, E)$, so we utilize Gauss' law for our final equation of motion

\nabla \cdot \mathbf{E} = \frac{\rho_q}{\epsilon_0} \implies \frac{\partial E}{\partial x} = \frac{e}{\epsilon_0}\left(n_j - n_e\right)

where $\rho_q$ is the charge density of the material.

Linearizing the Equations

We first find an equilibrium around which we can expand. Let the equilibrium state of our system be represented by $(n_{e0}, v_{e0} n_{i0}, v_{i0}, E_0)$. We assume that our equillibrium is steady in steady in state $(\partial_t = 0)$. This assumption immediately implies $E_0 = 0$. Furthermore, we can apply it to our other equations and attain the following relations

\frac{\partial}{\partial x}(n_{i0} v_{i0}) = 0

v_{i0} \frac{\partial v_{i0}}{\partial x} = 0

\frac{e}{\epsilon_0}(n_{j0} - n_{e0}) = 0 \implies \boxed{n_{j0} = n_{e0} = n_0},

where $i$ represents either an electron or ion equation. From the first two relations we can see that we are free to choose our conditions for the zeroth order velocities. We choose the initial conditions to be these velocities, namely $v_{e0} = V_0$ and $v_{j0} = 0$.

We now have an equilibrium state, namely $(n_{e0}, v_{e0}, n_{i0}, v_{i0}, E_0) = \left(n_0, V_0, n_0, 0, 0\right)$, which we can expand around.

We express each state function as an expansion about our equilibrium in orders of $\delta$. This process results in the following relations.

n_e = n_0 + \delta n_{e1} + O(\delta^2)

n_j = n_0 + \delta n_{j1} + O(\delta^2)

v_e = V_0 + \delta v_{e1} + O(\delta^2)

v_j = \delta v_{j1} + O(\delta^2)

E = \delta E_{1} + O(\delta^2).

We disregard the order 2 and higher terms, and substitute what remains into our original set of 5 equations. During this substitution we also disregard terms of order 2 or higher, as seen below.

\frac{\partial n_{e1}}{\partial t} + n_0 \frac{\partial v_{e1}}{\partial x} + V_0 \frac{\partial n_{e1}}{\partial x} = 0

m_e n_0 \left( \frac{\partial v_{e1}}{\partial t} + V_0 \frac{\partial v_{e1}}{\partial x} \right) = - e n_0 E_1

\frac{\partial n_{j1}}{\partial t} + n_0 \frac{\partial v_{j1}}{\partial x} = 0

m_i n_0 \frac{\partial v_{j1}}{\partial t} = e n_0 E_1

\frac{\partial E_1}{\partial x} = \frac{e}{\varepsilon_0}\,(n_{j1} - n_{e1})

Solving for the Dispersion Relation

We now make the following ansatz: each first order term (say $u_1$) can be expressed as a plane wave, where

$$u_1 = \tilde{u}_1 e^{i\left(kx - \omega t\right)}, $$

where $\tilde{u}_1 \in \mathbb{C}$ is the amplitude, $k \in \mathbb{R}^+$ is the wavenumber, and $\omega \in \mathbb{C}$ is the angular frequency with imaginary component $\operatorname{Im}(\omega) = \gamma$.

This transformation sends $\partial_t \to - i \omega$ and $\partial_x \to ik$, which in turn converts our 5 linear equations into 5 linear algebraic equations, which are expressed below.

-i \omega \tilde{n}_{e1} + n_0 (i k \tilde{v}_{e1}) + V_0 (i k \tilde{n}_{e1}) = 0

m_e n_0 \left\{ -i \omega \tilde{v}_{e1} + V_0 (i k \tilde{v}_{e1}) \right\} = - e n_0 \tilde{E}_1

-i \omega \tilde{n}_{j1} + n_0 (i k \tilde{v}_{j1}) = 0

m_j n_0 (-i \omega \tilde{v}_{j1}) = e n_0 \tilde{E}_1

i k \tilde{E}_1 = \frac{e}{\varepsilon_0}(\tilde{n}_{j1} - \tilde{n}_{e1})

Now, to solve this set of equations we first represent them in standard form, where

\tilde{n}_{e1}\left(V_0 k- \omega\right)+\tilde{v}_{e1}\left(n_0 k\right)=0

\tilde{v}_{e1}\left(i m_e (V_0 k - \omega) \right)+\tilde{E}_1\left(e \right)=0

\tilde{n}_{j1}(- \omega)+\tilde{v}_{j1}\left(n_0 k\right)=0

\tilde{v}_{j1}\left(m_j i \omega\right)+\tilde{E}_1\left(e \right)=0

\tilde{n}_{e1}(1)+\tilde{n}_{j1}(-1)+\tilde{E}_1\left(\frac{\varepsilon_0}{e} i k\right)=0

Now, setting $\mathbf{u} = [\tilde{n}_{e1}, \tilde{n}_{j1}, \tilde{v}_{e1}, \tilde{v}_{j1}, \tilde{E}_1]$ allows us to write all our equations as a matrix product, where

.latex

\begin{pmatrix}
V_0  k -  \omega & 0 & n_0 k & 0 & 0 \\ 
0 & 0 & i m_e(V_0 k - \omega) & 0 & e \\ 
0 & - \omega & 0 & n_0 k & 0 \\ 
0 & 0 & 0 & i m_j \omega & e  \\ 
1 & -1 & 0 & 0 & (\epsilon_0 i k) / e 
\end{pmatrix} 
\begin{pmatrix}
\tilde{n}_{e1} \\ \tilde{n}_{j1} \\ \tilde{v}_{e1} \\ \tilde{v}_{j1} \\ \tilde{E}_1
\end{pmatrix} = 0 \quad \iff \quad M \mathbf{u} = 0

This system has a unique solution for $\det(M) = 0$, which is the condition that gives us our dispersion relation. We implement this linear system in Mathematica and solve it below.

(*Define our variables*)
\[Epsilon]0 = Subscript[\[Epsilon], 0]; me = Subscript[m, e]; mj = Subscript[m, j]; n0 = Subscript[n, 0];
V0 = Subscript[V, 0];

(*Define our matrix M*)
M = {
    {V0 k -\[Omega], 0, k n0, 0, 0}, 
    {0, 0, I(V0 k - \[Omega])me, 0, e}, 
    {0, -\[Omega], 0, k n0, 0}, 
    {0, 0, 0, I mj \[Omega], e},
    {1, -1, 0, 0, (I \[Epsilon]0 k)/e}
}; MatrixForm[M]

(*Define and reduce dispersion equation*)
FullSimplify[Det[M] == 0]

From here we define our plasma frequencies and use them to simplify further.

\omega_{pe}^2 = \frac{n_0 e^2 }{\epsilon_0 m_e}, \quad \omega_{pj}^2 = \frac{n_0 e^2 }{\epsilon_0 m_j}

Notice that in the expression obtained by mathematica, we can cancel the $k$ values and multiply by $e$ to find

e^2 \omega^2 m_j n_0 + m_e \left(\omega - k V_0\right)^2 \left(e^2 n_0 - \omega^2 m_j \epsilon_0\right) = 0 .

Now from the plasma frequencies we have that $n_0 e^2 = \omega_{pe}^2 \epsilon_0 m_e = \omega_{pj}^2 \epsilon_0 m_j$, so we utilize this to write

\left(\omega_{pj}^2 \epsilon_0 m_j\right) \omega^2 m_j + m_e \left(\omega - k V_0\right)^2 \left(\omega_{pj}^2 \epsilon_0 m_j - \omega^2 m_j \epsilon_0\right) = 0.

We factor out an $m_j \epsilon_0$ term, and find that

m_j \omega_{pj}^2 \omega^2 + m_e\left(\omega - k V_0\right)^2 \left(\omega_{pj}^2 - \omega^2\right) = 0.

Now we divide both sides of the equation by $m_e$, and note that $m_j / m_e = \omega_{pe}^2 / \omega_{pj}^2$, so

\omega_{pe}^2 \omega^2 + \left(\omega - k V_0\right)^2 \left(\omega_{pj}^2 - \omega^2\right) = 0.

To make the structure clearer, we rearrange the equation by moving the second term to the right side and absorbing the negative sign.

\omega_{pe}^2 \omega^2 = \left(\omega - k V_0\right)^2 \left(\omega^2 - \omega_{pj}^2\right)

To normalize, we divide the entire equation by $\omega_{pe}^4$, which will allow us to form the dimensionless groups.

\frac{\omega_{pe}^2 \omega^2}{\omega_{pe}^4} = \frac{\left(\omega - k V_0\right)^2}{\omega_{pe}^2} \frac{\left(\omega^2 - \omega_{pj}^2\right)}{\omega_{pe}^2}

Now we can group the terms to explicitly match the normalized definitions.

\left(\frac{\omega}{\omega_{pe}}\right)^2 = \left(\frac{\omega}{\omega_{pe}} - \frac{k V_0}{\omega_{pe}}\right)^2 \left(\frac{\omega^2}{\omega_{pe}^2} - \frac{\omega_{pj}^2}{\omega_{pe}^2}\right)

Substituting the definitions $\Omega = \omega / \omega_{pe}$, $X = k V_0 / \omega_{pe}$, and $\mu^2 = \omega_{pj}^2 / \omega_{pe}^2$, we get the compact form.

\Omega^2 = (\Omega - X)^2 (\Omega^2 - \mu^2)

Finally, we arrange this into the standard polynomial form used for our computational analysis.

(\Omega^2 - \mu^2)(\Omega - X)^2 - \Omega^2 = 0.

(* Define new variables and our equation *)
eq2 = (\[CapitalOmega]^2 - \[Mu]^2)(\[CapitalOmega] - X)^2 - \[CapitalOmega]^2 == 0;

roots[x_] := \[CapitalOmega]/. Solve[eq2/.{\[Mu] -> 1/(Sqrt[1836]), X -> x}, \[CapitalOmega]]

Plot[{Im[roots[x][[1]]], Im[roots[x][[2]]],
          Im[roots[x][[3]]],Im[roots[x][[4]]]}, {x, 0, 2},
          PlotLegends -> {"Root 1", "Root 2", "Root 3", "Root 4"},
          AxesLabel -> {"X","Im(\[CapitalOmega])"}, PlotRange -> All ]

It seems clear that the root 3 is the instability, (as root 2 is exponential decay) but for what range of $X$ is the imaginary value of $\Omega$ nonzero? Next we use a numerical function to find this range.

(* --- Define a "Growth Rate" Function --- *)
growthRate[x_?NumericQ] := Max[Im[roots[x]]]


(* --- Use FindRoot to Find Where the Growth Rate is Zero --- *)



Print["--- Precise Instability Boundary ---"];

criticalXResult = FindRoot[
  growthRate[x] == 0,  
  {x, 1.1}             (* Our starting guess for the value of x *)
]

x_crit = x /. criticalXResult;

{x->1.1573469675922023`}

Let this value be $x_c$. From above we see that the instability is valid only for $0 \leq X \leq x_c$. Expanding our definitions, we have that instabilities only form when

0 \leq \frac{k v_0}{\omega_{pe}} \leq x_c \implies k \leq \frac{x_c \omega_{pe}}{v_0}.

Given a simulation of size $L$, the only waves that can exist have integer multiples $L$ as their wavelength. That is,

L = \lambda n \implies L = \frac{2 \pi n}{k} \implies k = \frac{2 \pi n}{L}.

Substituting this into our previous condition, we can find the minimal length of a simulation in which instabilites will be present, where setting $n=1$ gives us

\frac{2\pi n}{L} \leq \frac{x_c \omega_{pe}}{v_0} \implies \boxed{L \geq \frac{2 \pi v_0}{x_c \omega_{pe}}}